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Efficient implementation of the 
Hardy- Ramanuj an-Rademacher formula 

Fredrik Johansson 

o 

Abstract 

We describe how the Hardy- Ramanuj an-Rademacher formula can be implemented to allow the 
partition function p(n) to be computed with softly optimal complexity 0(n 1//2+o( ^) and very 
little overhead. A new implementation based on these techniques achieves speedups in excess of 
a factor 500 over previously published software and has been used by the author to calculate 
p(10 19 ), an exponent twice as large as in previously reported computations. 

We also investigate performance for mult i- evaluation of p(n), where our implementation of the 
Hardy- Ramanuj an-Rademacher formula becomes superior to power series methods on far denser 
sets of indices than previous implementations. As an application, we determine over 22 billion 
, new congruences for the partition function, extending Weaver's tabulation of 76,065 congruences. 
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1. Introduction 



£\\ ' Let p(n) denote the number of partitions of n, or the number of ways that n can be written as 

a sum of positive integers without regard to the order of the terms (A000041 in [OEIS ). The 
classical way to compute p(n) uses the generating function representation of p(n) combined 
with Euler's pentagonal number theorem 



<3\ 

in 

in 
o _ 

from which one can construct the recursive relation 



= II rh = £ (-D fe * fe(3fc - 1)/2 (i.i) 

n—0 k—1 \k= — oo / 



p(n) = £(-l) fc+1 (p(r 



fc(3fc - 1) V P fn- fc(3fc + 1) )). (1.2) 



k=l 

Equation (jl .2 j) provides a simple and reasonably efficient way to compute the list of values 
p(0),p(l), . . . ,p(n — l),p(n). Alternatively, applying FFT-based power series inversion to the 
right-hand side of ([1.1)) gives an asymptotically faster, essentially optimal algorithm for the 
same set of values. 

An attractive feature of Euler's method, in both the recursive and FFT incarnations, is that 
the values can be computed more efficiently modulo a small prime number. This is useful for 
investigating partition function congruences, such as in a recent large-scale computation of 
p{n) modulo small primes for n up to 10 9 [CDJPS07]. 

While efficient for computing p(n) for all n up to some limit, Euler's formula is impractical 
for evaluating p(n) for an isolated, large n. One of the most astonishing number-theoretical 
discoveries of the 20th century is the Hardy- Ramanuj an-Rademacher formula, first given as an 
asymptotic expansion by Hardy and Ramanuj an in 1917 [HR18] and subsequently refined 
to an exact representation by Rademacher in 1936 [Ra38], which provides a direct and 
computationally efficient expression for the single value p{n). 
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Simplified to a first-order estimate, the Hardy- Ramanuj an- Rademacher formula states that 



p(n) 



4nV3 



(1.3) 



from which one gathers that p(n) is a number with roughly n 1 / 2 decimal digits. The full version 
can be stated as 



<>»U [<M)+R( n ,N), 



»<»>=x:(fe)- 

U(x) = cosh(x) - Smh ^ , C(n) = ^V24n-1, 

4W = 5^* g cd(/i,fe),iexp (ni 
h=o ^ 



/c) 
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where s(h : k) is the Dedekind sum 



i=l 



and where the remainder satisfies |i?(n, 7V)| < M(n,N) with 



44tt 2 „_ 1/9 . tt\/2 / AT ^ 1/2 
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(1.4) 
(1.5) 
(1.6) 

(1.7) 

(1.8) 



It is easily shown that M^cn 1 / 2 ) ~ n -1 / 4 for every positive c. Rademacher's bound (jl.8p 
therefore implies that 0(n x l 2 ) terms in ([1.4)) suffice to compute p(n) exactly by forcing 
|i?(n, N)\ < 1/2 and rounding to the nearest integer. For example, we can take N = [n 1 / 2 ] 
when n > 65. 

In fact, it was pointed out by Odlyzko |Od95[ lKn05 that the Hardy- Ramanuj an- 
Rademacher formula "gives an algorithm for calculating p(n) that is close to optimal, since 
the number of bit operations is not much larger than the number of bits of p(n)". In other 
words, the time complexity should not be much higher than the trivial lower bound ^(n 1 / 2 ) 
derived from ([1.3)) just for writing down the result. Odlyzko's claim warrants some elaboration, 
since the Hardy- Ramanuj an- Rademacher formula ostensibly is a triply nested sum containing 
0(n 3 / 2 ) inner terms. 

The computational utility of the Hardy- Ramanuj an- Rademacher formula was, of course, 
realized long before the availability of electronic computers. For instance, Lehmer [Le36j used 
it to verify Ramanujan's conjectures p(599) = mod 5 3 and £>(721) = mod ll 2 . Implementa- 
tions are now available in numerous mathematical software systems, including Pari/GP [Pari], 
Mathematica [Woll] and Sage Sage|. However, apart from Odlyzko's remark, we find few 
algorithmic accounts of the Hardy- Ramanuj an- Rademacher formula in the literature, nor any 
investigation into the optimality of the available implementations. 

The present paper describes a new C implementation of the Hardy- Ramanuj an- Rademacher 
formula. The code is freely available as a component of the Fast Library for Number Theory 
(FLINT) [HalO], released under the terms of the GNU General Public License. We show 
that the complexity for computing p(n) indeed can be bounded by 0(n 1 / 2+o( ^ 1 ^), and observe 
that our implementation comes close to being optimal in practice, improving on the speed of 
previously published software by more than two orders of magnitude. 

We benchmark the code by computing some extremely large isolated values of p(n). We also 
investigate efficiency compared to power series methods for evaluation of multiple values, and 
finally apply our implementation to the problem of computing congruences for p(n). 
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2. Simplification of exponential sums 

A naive implementation of formulas ()1 .4|) — ()1 .7|) requires 0(n 3 / 2 ) integer operations to 
evaluate Dedekind sums, and 0(n) numerical evaluations of complex exponentials (or cosines, 
since the imaginary parts ultimately cancel out). In the following section, we describe how the 
number of integer operations and cosine evaluations can be reduced, for the moment ignoring 
numerical evaluation. 

A first improvement, used for instance in the Sage implementation, is to recognize that 
Dedekind sums can be evaluated in 0(\ogk) steps using a GCD-style algorithm, as described 



by Apostol | Ap97 , or with Knuth's fraction- free algorithm [Kn77 which avoids the overhead 
of rational arithmetic. This reduces the total number of integer operations to 0(n log n), which 
is a dramatic improvement but still leaves the cost of computing p(n) quadratic in the size of 
the final result. 

Fortunately, the Ak(n) sums have additional structure as discussed in [Le37, Le38l |RW41| 
Wh56, Ha70 , allowing the computational complexity to be reduced. Since numerous imple- 
menters of the Hardy- Ramanuj an- Rademacher formula until now appear to have overlooked 
these results, it seems appropriate that we reproduce the main formulas and assess the 
computational issues in more detail. We describe two concrete algorithms: one simple, and 
one asymptotically fast, the latter being implemented in FLINT. 



2.1. A simple algorithm 

Using properties of the Dedekind eta function, one can derive the formula (which Whiteman 
|Wh56j attributes to Selberg) 



A k (n) = 




in which the summation ranges over < I < 2k and only 0(k 1 ^ 2 ) terms are nonzero. With a 
simple brute force search for solutions of the quadratic equation, this representation provides 
a way to compute Ak(n) that is both simpler and more efficient than the usual definition (|1.6|) . 

Although a brute force search requires 0(k) loop iterations, the successive quadratic terms 
can be generated without multiplications or divisions using two coupled linear recurrences. 
This only costs a few processor cycles per loop iteration, which is a substantial improvement 
over computing Dedekind sums, and means that the cost up to fairly large k effectively will 
be dominated by evaluating 0(k 1 ^ 2 ) cosines, adding up to 0(n 3 / 4 ) function evaluations for 
computing p(n). 

A basic implementation of (j2.ip is given as Algorithm [TJ Here the variable m runs over 
the successive values of (3/ 2 +Z)/2, and r runs over the differences between consecutive m. 
Various improvements are possible: a modification of the equation allows cutting the loop 
range in half when k is odd, and the number of cosine evaluations can be reduced by counting 
the multiplicities of unique angles after reduction to [0,7r/4), evaluating a weighted sum 
J2wicos(0i) at the end - possibly using trigonometric addition theorems to exploit the fact 
that the differences 6i+i — Oi between successive angles tend to repeat for many different i. 



2.2. A fast algorithm 

From Selberg's formula (|2.ip . a still more efficient but considerably more complicated 
multiplicative decomposition of Ak(n) can be obtained. The advantage of this representation 
is that it only contains O(logfc) cosine factors, bringing the total number of cosine evaluations 
for p(n) down to 0(n x / 2 logn). It also reveals exactly when Ak{n) = (which is about half the 
time). We stress that these results are not new; the formulas are given in full detail and with 
proofs in [Wh56]. 
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Algorithm 1 Simple algorithm for evaluating Ak(n) 
Input: Integers fc, n > 

Output: s = Ak(n), where Ak(n) is defined as in (|1.6|) 
if fc < 1 then 

return fc 
else if = 2 then 

return (— l) n 
end if 

(s, r, m) ^— (0, 2, (n mod fc)) 
for < Z < 2fc do 
if m = then 

s + (-1)' cos (tt(6Z + l)/(6fc)) 

end if 

m ^— m + r 

if m > fc then m ^— m — k {m ^— rn mod fc} 

r «— r + 3 

if r > fc then r r — fc {r ^— r mod fc} 
end for 

return (fc/3) 1 / 2 5 



First consider the case when fc is a power of a prime. Clearly A\(n) — 1 and A<i(n) — (— l) n . 
Otherwise let fc = p x and v = 1 — 24n. Then, using the notation (a|m) for Jacobi symbols to 
avoid confusion with fractions, we have 



'(-l) A (-l|ra 2 )fc 1/2 sin(47rra 2 /8fc) 
A k {n) = I 2(-l) A + 1 (m 3 |3)(fc/3) 1 / 2 sin(47rm 3 /3fc) 
^2(3|fc)fc 1 / 2 cos(47rm p /fc) 

where m 2 , ms and m p respectively are any solutions of 

(3ra 2 ) 2 = v mod 8fc 
(Sms) 2 = v mod 3fc 
(24m p ) 2 = v mod fc 



ifp = 2 
if p = 3 
ifp>3 



(2.2) 



(2.3) 
(2.4) 
(2.5) 



provided, when p > 3, that such an m p exists and that gcd(v, fc) = 1. If, on the other hand, 
p > 3 and either of these two conditions do not hold, we have 



A k (ri) = < 





(3|fc)fc x / 2 




if v is not a quadratic residue modulo fc 
if v = mod p, A = 
if v = mod p, A > 1. 



(2.6) 



If fc is not a prime power, assume that fc = fcifc 2 where gcd(fci, fc 2 ) = 1. Then we can factor 
Ak(n) as Ak(n) = A kl (ni)Ak 2 (^2), where ni,n 2 are any solutions of the following equations. 
If fci = 2, then 



32n 2 



8n + 1 mod fc 2 

*2 



ni =n-(fc|-l)/8 mod 2, 



if fci = 4, then 



Jl28n 2 = 8n + 5 mod fc 2 

|fc|ni = n - 2 - (fc| - l)/8 mod 4, 



(2.7) 



(2.8) 
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and if k\ is odd or divisible by 8, then 

/cfc^eni = d 2 en + (k 2 — l)/^i mod k\ ^ ^ 

k\d\en 2 = (iien + (/c 2 — l)/c?2 mod &2 

where di = gcd(24, hi), d 2 = gcd(24, fe 2 ), 24 = did 2 e. 

Here (/c 2 — l)/d denotes an operation done on integers, rather than a modular division. All 
other solving steps in ([2.2p ~ (j2.9p amount to computing greatest common divisors, carrying 
out modular ring operations, finding modular inverses, and computing modular square roots. 
Repeated application of these formulas results in Algorithm [2j where we omit the detailed 
arithmetic for brevity. 



Algorithm 2 Fast algorithm for evaluating Akin) 
Input: Integers k > 1, n > 

Output: s = Ak(n), where Ak(n) is defined as in (|1.6|) 
Compute the prime factorization k = p^p 2 2 • • - Pj 3 
s <- 1 

for 1 < i < j and while s ^ do 
if i < j then 

Compute ni,n2 by solving the respective case of (|2.7p - ([2.9|) 
s ^— s x Ak 1 (ni) {Handle the prime power case using (j2.2p ~ (|2.6p } 
(fe,n) <- (k 2 ,n 2 ) 
else 

s <— s x Ak{n) {Prime power case} 
end if 
end for 
return «s 



2.3. Computational cost 

A precise complexity analysis of Algorithm [2] should take into account the cost of integer 
arithmetic. Multiplication, division, computation of modular inverses, greatest common divisors 
and Jacobi symbols of integers bounded in absolute value by 0(k) can all be performed with 
bit complexity 0(log 1+o(1) k). 

At first sight, integer factorization might seem to pose a problem. We can, however, factor 
all indices k summed over in (|1.4|) in 0(n x / 2 log 1+o( ^ n) bit operations. For example, using 
the sieve of Eratosthenes, we can precompute a list of length n 1 / 2 where entry k is the largest 
prime dividing k. 

A fixed index k is a product of at most O(logfc) prime powers with exponents bounded 
by O(logfc). For each prime power, we need 0(1) operations with roughly the cost of 
multiplication, and 0(1) square roots, which are the most expensive operations. 

To compute square roots modulo p A , we can use the Tonelli- Shanks algorithm [T o91[[Sh72j 
or Cipolla's algorithm |Ci03j modulo p followed by Hensel lifting up to p x . Assuming that 
we know a quadratic nonresidue modulo p, the Tonelli- Shanks algorithm requires 0(log 3 k) 
multiplications in the worst case and 0(log 2 k) multiplications on average, while Cipolla's 
algorithm requires 0(log 2 k) multiplications in the worst case [CP05 . This puts the bit 
complexity of factoring a single exponential sum A k {n) at O(log 3+o(1) jfe), and gives us the 
following result: 
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Theorem 1 . Assume that we know a quadratic nonresidue modulo p for all primes p up to 
n I . Then we can factor all the Ak(n) required for evaluating p(n) using 
bit operations. 

The assumption in Theorem [I] can be satisfied with a precomputation that does not affect 
the complexity. If n2(pk) denotes the least quadratic nonresidue modulo the kth prime number, 
it is a theorem of Erdos [Er6lJ IPol2j that as x — >• oo, 

oo 

E n M "> E | = C < 3 ' 675 - ( 2 - 10 ) 

Pk<x k=l 

Given the primes up to x = n 1 / 2 , we can therefore build a table of nonresidues by testing 
no more than (C + o(l))7r(n 1 / 2 ) candidates. Since 7r(n 1//2 ) = 0{in}^ 2 / \ogn) and a quadratic 
residue test takes 0(log 1+0( ^ p) time, the total precomputation time is Oin 1 ! 2 log 0( ^ n). 

In practice, it is sufficient to generate nonresidues on the fly since 0(1) candidates need to 
be tested on average, but we can only prove an 0(log c k) bound for factoring an isolated Ak(n) 
by assuming the Extended Riemann Hypothesis which gives n2(p) — 0(log 2 p) [An52 . 

2.4. Implementation notes 

As a matter of practical efficiency, the modular arithmetic should be done with as little 
overhead as possible. FLINT provides optimized routines for arithmetic with moduli smaller 
than 32 or 64 bits (depending on the hardware word size) which are used throughout; including, 
among other things, a binary-style GCD algorithm, division and remainder using precomputed 
inverses, and supplementary code for operations on two-limb (64 or 128 bit) integers. 

We note that since Ak(n) = Ak(n + fc), we can always reduce n modulo fc, and perform 
all modular arithmetic with moduli up to some small multiple of k. In principle, the 
implementation of the modular arithmetic in FLINT thus allows calculating p(n) up to 
approximately n = (2 64 ) 2 w 10 38 on a 64-bit system, which roughly equals the limit on n 
imposed by the availability of addressable memory to store p(n). 

At present, our implementation of Algorithm [2] simply calls the FLINT routine for integer 
factorization repeatedly rather than sieving over the indices. Although convenient, this 
technically results in a higher total complexity than 0(n 1 / 2+o( ^ 1 ^). However, the code for 
factoring single-word integers, which uses various optimizations for small factors and Hart's 
"One Line Factor" variant of Lehman's method to find large factors [Hall , is fast enough 
that integer factorization only accounts for a small fraction of the running time for any feasible 
n. If needed, full sieving could easily be added in the future. 

Likewise, the square root function in FLINT uses the Tonelli- Shanks algorithm and generates 
a nonresidue modulo p on each call. This is suboptimal in theory but efficient enough in practice. 



3. Numerical evaluation 

We now turn to the problem of numerically evaluating (|1.4|) - (|1.5|) using arbitrary-precision 
arithmetic, given access to Algorithm[2]for symbolically decomposing the Ak(n) sums. Although 
(|1.8p bounds the truncation error in the Hardy- Ramanuj an- Rademacher series, we must also 
account for the effects of having to work with finite-precision approximations of the terms. 

3.1. Floating-point precision 

We assume the use of variable-precision binary floating-point arithmetic (a simpler but less 
efficient alternative, avoiding the need for detailed manual error bounds, would be to use 
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arbitrary-precision interval arithmetic). Basic notions about floating-point arithmetic and error 
analysis can be found in [Hi02 . 

If the precision is r bits, we let e — 2~ r denote the unit roundoff. We use the symbol x 
to signify a floating-point approximation of an exact quantity x, having some relative error 
S = (x — x)/x when x ^ 0. If x is obtained by rounding x to the nearest representable floating- 
point number (at most 0.5 ulp error) at precision r, we have \5\ < e. Except where otherwise 
noted, we assume correct rounding to nearest. 

A simple strategy for computing p(n) is as follows. For a given n, we first determine an 
N such that \R(n,N)\ < 0.25, for example using a linear search. A tight upper bound for 
log 2 M{n,N) can be computed easily using low-precision arithmetic. We then approximate the 
fcth term tk using a working precision high enough to guarantee 

\tk~t k \ < — , (3.1) 

and perform the outer summation such that the absolute error of each addition is bounded by 
0.125/7V. This clearly guarantees \p(n) — p(n)\ < 0.5, allowing us to determine the correct value 
of p{n) by rounding to the nearest integer. We might, alternatively, carry out the additions 
exactly and save one bit of precision for the terms. 

In what follows, we derive a simple but essentially asymptotically tight expression for a 
working precision, varying with /c, sufficiently high for (|3.1|) to hold. Using Algorithm [2j we 
write the term to be evaluated in terms of exact integer parameters a, /3, a, b,pi, qi as 




Lemma 2. Letp G Z, q G N + and let r be a precision in bits with 2 r > max(3g, 64). Suppose 
that sin and cos can be evaluated on (0,7r/4) with relative error at most 2s for floating-point 
input, and suppose that tt can be approximated with relative error at most e. Then we can 
evaluate cos(p7r/q) with relative error less than 5. be. 



Proof. 

We first reduce p and q with exact integer operations so that < Ap < q, giving an angle 
in the interval (0,7r/4). Then we approximate x = pir/q using three roundings, giving x = 
x(l + 5 X ) where |^| < (1 + e) 3 - 1. The assumption e < l/(Sq) gives (q/(q - 1))(1 + 5 X ) < 1 
and therefore also x G (0,7r/4). 

Next, we evaluate f(x) where / = ± cos or / = zb sin depending on the argument reduction. 
By Taylor's theorem, we have f(x)=f(x)(l + 5' x ) where 

,,,, \m-f(x)\ x\s x \\f(o\ 

lsJ = m = (3 - 3) 

for some £ between x and x : giving \S f x \ < (j7ry/2)\5 x \. Finally, rounding results in 

fix) = f(x)(l + S) = f(x)(l + S' x )(l + 5 f ) 
where \Sf\ < 2e. The inequality e < 1/64 gives \S\ < b.be. □ 

To obtain a simple error bound for U(x) where x = C/k, we make the somewhat crude 
restriction that n > 2000. We also assume k < n 1 / 2 and x > 3, which are not restrictions: if N 
is chosen optimally using Rademacher's remainder bound (|1.8)h the maximum k decreases and 
the minimum x increases with larger n. In particular, n > 2000 is sufficient with Rademacher's 
bound (or any tighter bound for the remainder). 
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We assume that C is precomputed; of course, this only needs to be done once during the 
calculation of p(n), at a precision a few bits higher than that of the k = 1 term. 



Lemma 3. Suppose n > 2000 and let r be a precision in bits such that 2 r > 
max(16n 1 / 2 , 2 10 ). Let x = C/k where C is defined as in (|1.5|) and where k is constrained such 
that k < n 1 / 2 and x > 3. Assume that C = C(n)(\ + 8 c) has been precomputed with \Sc\ < 2e 
and that sinh and cosh can be evaluated with relative error at most 2s for floating-point input. 
Then we can evaluate U(x) with relative error at most (9x + 15) e. 



Proof We first compute x = x(l + 8 X ) — (C/k)(l + 5c)(l + #o) where \So\ < e. Next, we 
compute 

U(x) = U(x)(l + Su) = U(x)(l + 0(1 + fo) = U(x)(l + J) (3.4) 

where we have to bound the error 8' x propagated in the composition as well as the rounding 
error 8u in the evaluation of U(x). Using the inequality x\S x \ < Axe < log 2, we have 

iw i / x\6 x \U'(x + x\6 x \) x\6 x \exp(x + x\6 x \) x\5 x \exp(x) . , 

^ - ^) - 2t^) - - (3 - 5) 

Evaluating U (x) using the obvious sequence of operations results in 



\&u\ = 



cosh(£)(l + 2Ji) - + 2J 2 )(1 + S 3 ) ) (1 + J 4 ) - U(x) 



(3.6) 



/7(f) 

where \8i\ < e and x > z where z = 3(1 — 4e). This expression is maximized by setting £ as 
small as possible and taking Si = 84 = — S2 = —£3 = e, which gives 

. r . cosh(z) , n . sinh(z) , n 9n 

1 C/| < ~t7(I) £( + ^ + 7u{z) € ^ + £ " ^ < ( ) 

Expanding (|3^4|) using ([33]) and (pTTj) gives |5| < £(5.5 + 9x + 56xe + 33xe 2 ). Finally, we 
obtain 5.5 + 56xe + 33xe 2 < 15 by a direct application of the assumptions. □ 

Put together, assuming floating-point implementations of standard transcendental functions 
with at most 1 ulp error (implying a relative error of at most 2e), correctly rounded arithmetic 
and the constant 7r, we have: 



Theorem 4. Let n > 2000. For (|3.ip to hold, it is sufficient to evaluate (13. 2p using a 
precision of r = max(log 2 N + log 2 |tfc| + log 2 (10x + 7m + 22) + 3, | log 2 n + 5, 11) bits. 



Proof. We can satisfy the assumptions of lemmas [2] and [3] In particular, 3q < 24k < 
24n x / 2 < 2 r . The top-level arithmetic operations in (|3.2|h including the square roots, amount to 
a maximum of m + 6 roundings. Lemmas [2] and [3] and elementary inequalities give the relative 
error bound 

|5| < (1 + e) m+6 (1 + 5.5e) m (1 + (15 + 9x)e) - 1 (3.8) 

< (1 + 1 (r7 ; +6) ^ ( 1 + t^-) a + ( is + w - 1 ( 3 - 9 ) 

V l-(ra + 6)e/V l-5.5mey v v yy v ; 
_ 21s + 6.5m£ — 33ra£ 2 — 5.5ra 2 £ 2 + 9x£ 

(l-5.5era)(l-e(ra + 6)) ( ' ^ 

< (lOx + 7m + 22)e. (3.11) 
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The result follows by taking logarithms in (|3.ip . □ 

To make Theorem [H effective, we can use m < log 2 k and bound \tk \ using (|1.4p with < fc 
and ?7(x) < e x /2, giving 

, , , (24n-l) 1 / 2 7r logfc ^ _ / o log3\ 

log < i ^ + -f- - log(24n - 1) + (\og2 + -|-J . (3.12) 

Naturally, for n < 2000, the same precision bound can be verified to be sufficient through 
direct computation. We can even reduce overhead for small n by using a tighter precision, 
say r = \tk\ +0(1), up to some limit small enough to be tested exhaustively (perhaps much 
larger than 2000). The requirement that r > \ log 2 n + 0(1) always holds in practice if we set 
a minimum precision; for n feasible on present hardware, it is sufficient to never drop below 
IEEE double (53-bit) precision. 



3.2. Computational cost 

We assume that r-bit floating-point numbers can be multiplied in time M(r) — 
0(r log 1+0( ^ r). It is well known (see [BZllj ) that the elementary functions exp, log, sin 
etc. can be evaluated in time 0(M(r) logr) using methods based on the arithmetic-geometric 
mean (AGM). A popular alternative is binary splitting, which typically has cost 0(M(r) log 2 r) 
but tends to be faster than the AGM in practice. 

To evaluate p(n) using the Hardy-Ramanujan-Rademacher formula, we must add 0(n 1//2 ) 
terms each of which can be written as a product of O(logfc) factors. According to (j3.12ft and 
the error analysis in the previous section, the kth term needs to be evaluated to a precision of 
0(n x / 2 //c) + O(logn) bits. Using any combination of 0(M(r) log" r) algorithms for elementary 
functions, the complexity of the numerical operations is 

( nl//2 / l/ 2 \ 1/2 \ 

J2 log k M lo g" V ) = ° ( nV2 lo S a+3+ ° (1) n ) ( 3 ' 13 ) 

which is nearly optimal in the size of the output. Combined with the cost of the factoring 
stage, the complexity for the computation of p(n) as a whole is therefore, when properly 
implemented, softly optimal at 0(n 1 / 2+o( ^ 1 ^). From (j3. 13j) with the best known complexity 
bound for elementary functions, we obtain: 

Theorem 5. The value p(n) can be computed in time 0(n x / 2 log 4+0( ^ n). 

A subtle but crucial detail in this analysis is that the additions in the main sum must be 
implemented in such a way that they have cost 0(n 1 / 2 //c) rather than 0(n x / 2 ), since the latter 
would result in an 0{n) total complexity. If the additions are performed in-place in memory, 
we can perform summations the natural way and rely on carry propagation terminating in an 
expected 0(1) steps, but many implementations of arbitrary-precision floating-point arithmetic 
do not provide this optimization. 

One way to solve this problem is to add the terms in reverse order, using a precision that 
matches the magnitude of the partial sums. Or, if we add the terms in forward order, we can 
amortize the cost by keeping separate summation variables for the partial sums of terms not 
exceeding 7*1, r*i/2, r*i/4, 7*1/8, . . . bits. 
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3.3. Arithmetic implementation 

FLINT uses the MPIR library, derived from GMP, for arbitrary-precision arithmetic, and 
the MPFR library on top of MPIR for asymptotically fast arbitrary-precision floating-point 
numbers and correctly rounded transcendental functions [MPIR, GMP , MPFR . Thanks to 
the strong correctness guarantees of MPFR, it is relatively straightforward to write a provably 
correct implementation of the partition function using Theorem |4j 

Although the default functions provided by MPFR are quite fast, order-of-magnitude 
speedups were found possible with custom routines for parts of the numerical evaluation. An 
unfortunate consequence is that our implementation currently relies on routines that, although 
heuristically sound, have not yet been proved correct, and perhaps are more likely to contain 
implementation bugs than the well-tested standard functions in MPFR. 

All such heuristic parts of the code are, however, well isolated, and we expect that they can 
be replaced with rigorous versions with equivalent or better performance in the future. 

3.4. Hardware arithmetic 

Inspired by the Sage implementation, which was written by Jonathan Bober, our implemen- 
tation switches to hardware (IEEE double) floating-point arithmetic to evaluate (|3.2p when 
the precision bound falls below 53 bits. This speeds up evaluation of the "long tail" of terms 
with very small magnitude. 

Using hardware arithmetic entails some risk. Although the IEEE floating-point standard 
implemented on all modern hardware guarantees 0.5 ulp error for arithmetic operations, 
accuracy may be lost, for example, if the compiler generates long-double instructions which 
trigger double rounding, or if the rounding mode of the processor has been changed. 

We need to be particularly concerned about the accuracy of transcendental functions. 
The hardware transcendental functions on the Intel Pentium processor and its descendants 
guarantee an error of at most 1 ulp when rounding to nearest [In95 , as do the software 
routines in the portable and widely used FDLIBM library [Sun]. Nevertheless, some systems 
may be equipped with poorer implementations. 

Fortunately, the bound (|1.8|) and Theorem [4] are lax enough in practice that errors up 
to a few ulp can be tolerated, and we expect any reasonably implemented double-precision 
transcendental functions to be adequate. Most importantly, range reducing the arguments of 
trigonometric functions to (0,7r/4) avoids catastrophic error for large arguments which is a 
misfeature of some implementations. 

3.5. High-precision evaluation of exponentials 

MPFR implements the exponential and hyperbolic functions using binary splitting at high 
precision, which is asymptotically fast up to logarithmic factors. We can, however, improve 
performance by not computing the hyperbolic functions in U(x) from scratch when k is 
small. Instead, we precompute exp(C) with the initial precision of C, and then compute 
(cosh(C/fc), sinh(C/fc)) from (exp(C)) 1 / /c ; that is, by fcth root extractions which have cost 
0((log k)M{r)). Using the builtin MPFR functions, root extraction was found experimentally 
to be faster than evaluating the exponential function up to approximately k = 35 over a large 
range of precisions. 

For extremely large n, we also speed up computation of the constant C by using binary 
splitting to compute tt (adapting code written by H. Xue |GMP2j ) instead of the default 
function in MPFR, which uses arithmetic-geometric mean iteration. As has been pointed out 
previously |Zi06 , binary splitting is more than four times faster for computing tt in practice, 
despite theoretically having a log factor worse complexity. When evaluating p(n) for multiple 
values of n, the value of tt should of course be cached, which MPFR does automatically. 
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3.6. High-precision cosines 

The MPFR cosine and sine functions implement binary splitting, with similar asymptotics 
as the exponential function. At high precision, our implementation switches to custom code 
for evaluating a = cos(p7r/q) when q is not too large, taking advantage of the fact that a is an 
algebraic number. Our strategy consists of generating a polynomial P such that P(a) = and 
solving this equation using Newton iteration, starting from a double precision approximation 
of the desired root. Using a precision that doubles with each step of the Newton iteration, the 
complexity is 0(deg(P)M(r)). 

The numbers cos(p7r/q) are computed from scratch as needed: caching values with small p and 
q was found to provide a negligible speedup while needlessly increasing memory consumption 
and code complexity. 

Our implementation uses the minimal polynomial $> n (x) of cos(27r/n), which has degree 
d = (j)(n)/2 for n > 3 [WZ93 . More precisely, we use the scaled polynomial 2 d <&{x) G Z[x]. This 
polynomial is looked up from a precomputed table when n is small, and otherwise is generated 
using a balanced product tree, starting from floating-point approximations of the conjugate 
roots. As a side remark, this turned out to be around a thousand times faster than computing 
the minimal polynomial with the standard commands in either Sage or Mathematica. 

We sketch the procedure for high-precision evaluation of cos(j>7r/g) as Algorithm [3j omitting 
various special cases and implementation details (for example, our implementation performs the 
polynomial multiplications over Z[x] by embedding the approximate coefficients as fixed-point 
numbers). 



Algorithm 3 High-precision numerical evaluation of cos(p7r/q) 
Input: Coprime integers p and q with q > 3, and a precision r 
Output: An approximation of cos(p7r/q) accurate to r bits 

n ^— (1 + (p mod 2)) q 

d <r- 0(n)/2 

{Bound coefficients in 2 d Yli=i( x ~ a )} 



Compute precisions ro = r + 8, r\ = ro/2 + 8, . . . , rj = rj-i/2 + 8 < 50 
x <— cos(p7r/q) {To 50 bits, using basecase algorithm} 
for i j — 1, j — 2 . . . do 

{Evaluate using the Horner scheme at n + b bit precision} 

x <- x - $>(x)/<fr'(x) 
end for 
return x 



We do not attempt to prove that the internal precision management of Algorithm [3] is correct. 
However, the polynomial generation can easily be tested up to an allowed bound for and 
the function can be tested to be correct for all pairs p, q at some fixed, high precision r. We 
may then argue heuristically that the well-behavedness of the object function in the root- 
finding stage combined with the highly conservative padding of the precision by several bits 
per iteration suffices to ensure full accuracy at the end of each step in the final loop, given an 
arbitrary r. 
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A different way to generate $> n (x) using Chebyshev polynomials is described in [WZ93 . 
One can also use the squarefree part of an offset Chebyshev polynomial 

Fix) - ^ {X) ~ 1 



gcd(T 2 ,(a;)-l,r^(a;)) 

directly, although this is somewhat less efficient than the minimal polynomial. 

Alternatively, since cos(p7r/q) is the real part of a root of unity, the polynomial x q — 1 could 
be used. The use of complex arithmetic adds overhead, but the method would be faster for large 
q since x q can be computed in time 0((\ogq)M(r)) using repeated squaring. We also note that 
the secant method could be used instead of the standard Newton iteration in Algorithm [3l This 
increases the number of iterations, but removes the derivative evaluation, possibly providing 
some speedup. 

In our implementation, Algorithm [3] was found to be faster than the MPFR trigonometric 
functions for q < 250 roughly when the precision exceeds 400 + Aq 2 bits. This estimate includes 
the cost of generating the minimal polynomial on the fly. 



3.7. The main algorithm 

Algorithm [4] outlines the main routine in FLINT with only minor simplifications. To avoid 
possible corner cases in the convergence of the HRR sum, and to avoid unnecessary overhead, 
values with n < 128 (exactly corresponding to p(n) < 2 32 ) are looked up from a table. We only 
use k, n, N in Theorem [4] in order to make the precision decrease uniformly, allowing amortized 
summation to be implemented in a simple way. 



Algorithm 4 Main routine implementing the HRR formula 
Input: n > 128 
Output: p(n) 

Determine N and initial precision r\ using Theorem [4] 

C <r- f V24n- 1 {At n + 3 bits} 

u <— exp(C) 

s\ <- s 2 <- 

for 1 < k < N do 

Write term k as (|3.2p by calling Algorithm [2] 
if A k {n) + then 

Determine term precision r& for \t^\ using Theorem [4] 

{Use Algorithm E] if q { < 250 and r k > 400 + 4<? 2 } 

t <- (-l) 8 y/a/bYlcos(pi7r/qi) 

t^tx U(C/k) {Compute U from u x l k if k < 35} 

{Amortized summation: r(s2) denotes precision of the variable S2} 

S2 <~ S 2 + t 

if 2r/e < r(s2) then 

$1 <— $1 + -§2 {Exactly or with precision exceeding n} 
^(^2) ^— {Change precision} 

S 2 ^0 

end if 
end if 
end for 

return [si + s 2 + 
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Since our implementation presently relies on some numerical heuristics (and in any case, 
considering the intricacy of the algorithm), care has been taken to test it extensively. All n < 10 6 
have been checked explicitly, and a large number of isolated n ^> 10 6 have been compared 
against known congruences and values computed with Sage and Mathematica. 

As a strong robustness check, we observe experimentally that the numerical error in the final 
sum decreases with larger n. For example, the error is consistently smaller than 10 -3 for n > 10 6 
and smaller than 10 -4 for n > 10 9 . This phenomenon reflects the fact that (|1.8|) overshoots 
the actual magnitude of the terms with large k, combined with the fact that rounding errors 
average out pseudorandomly rather than approaching worst-case bounds. 



4. Benchmarks 

Table [T] and Figure [T] compare performance of Mathematica 7, Sage 4.7 and FLINT on a 
laptop with a Pentium T4400 2.2 GHz CPU and 3 GiB of RAM, running 64-bit Linux. To the 
author's knowledge, Mathematica and Sage contain the fastest previously available partition 
function implementations by far. 

The FLINT code was run with MPIR version 2.4.0 and MPFR version 3.0.1. Since Sage 4.7 
uses an older version of MPIR and Mathematica is based on an older version of GMP, differences 
in performance of the underlying arithmetic slightly skew the comparison, but probably not 
by more than a factor two. 

The limited memory of the aforementioned laptop restricted the range of feasible n to 
approximately 10 16 . Using a system with an AMD Opteron 6174 processor and 256 GiB RAM 
allowed calculating p(10 17 ), p(10 18 ) and p(10 19 ) as well. The last computation took just less 
than 100 hours and used more than 150 GiB of memory, producing a result with over 11 billion 
bits. Some large values of p(n) are listed in Table [2j 

As can be seen in Tableland Figure [TJ the FLINT implementation exhibits a time complexity 
only slightly higher than 0(n 1 / 2 ), with a comparatively small constant factor. The Sage 



n 


Mathematica 7 


Sage 4.7 


FLINT 


Initial 


10 4 


69 ms 


1 ms 


0.20 ms 




10 5 


250 ms 


5.4 ms 


0.80 ms 




10 6 


590 ms 


41 ms 


2.74 ms 




10 7 


2.4 s 


0.38 s 


0.010 s 




10 8 


11 s 


3.8 s 


0.041 s 




10 9 


67 s 


42 s 


0.21 s 


43% 


10 10 


340 s 




0.88 s 


53% 


10 11 


2,116 s 




5.1 s 


48% 


IO 12 


10,660 s 




20 s 


49% 


IO 13 






88 s 


48% 


10 14 






448 s 


47% 


10 15 






2,024 s 


39% 


10 16 






6,941 s 


45% 


IO 17 






27,196* s 


33% 


10 18 






87,223* s 


38% 


IO 19 






350,172* s 


39% 



Table 1 . Timings for computing p(n) in Mathematica 7, Sage 4. 7 and FLINT up to 
n = 10 16 on the same system, as well as FLINT timings for n = 10 1T to 10 19 (*) done on 
different (slightly faster) hardware. Calculations running less than one second were repeated, 
allowing benefits from data caching. The rightmost column shows the amount of time in the 
FLINT implementation spent computing the first term. 
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10 5 F 1 1 1 1 1 1 1 1 1 1 r 




n 



Figure 1. CPU time t in seconds for computing p(n) : FLINT (blue squares), Mathematica 7 
(green circles), Sage 4.7 (red triangles). The dotted line shows t = lO -6 ?! 1 / 2 , indicating the 
slope of an idealized algorithm satisfying the trivial lower complexity bound ^(n 1 / 2 ) (the 

offset 10 -6 is arbitrary). 

implementation is fairly efficient for small n but has a complexity closer to O(n), and is limited 
to arguments n < 2 32 « 4 x 10 9 . 

The Mathematica implementation appears to have complexity slightly higher than 0(n 1 / 2 ) 
as well, but consistently runs about 200-500 times slower than our implementation. Based on 
extrapolation, computing p(10 19 ) would take several years. It is unclear whether Mathematica is 
actually using a nearly-optimal algorithm or whether the slow growth is just the manifestation 
of various overheads dwarfing the true asymptotic behavior. The ratio compared to FLINT 
appears too large to be explained by differences in performance of the underlying arithmetic 
alone; for example, evaluating the first term in the series for p(10 10 ) to required precision in 
Mathematica only takes about one second. 

We get one external benchmark from [BB10 , where it is reported that R. Crandall 
computed p(10 9 ) in three seconds on a laptop in December 2008, "using the Hardy- Ramanujan- 
Rademacher 'finite' series for p(n) along with FFT methods". Even accounting for possible 
hardware differences, this appears to be an order of magnitude slower than our implementation. 

4.1. Optimality relative to the first term 

Table [1] includes time percentages spent on evaluating the first term, exp(C), in the FLINT 
implementation. We find that this step fairly consistently amounts to just a little less than half 
of the running time. Our implementation is therefore nearly optimal in a practical sense, since 
the first term in the Hardy- Ramanuj an- Rademacher expansion hardly can be avoided and at 
most a factor two can be gained by improving the tail evaluation. 

Naturally, there is some potential to implement a faster version of the exponential function 
than the one provided by MPFR, reducing the cost of the first term. Improvements on the level 
of bignum multiplication would, on the other hand, presumably have a comparatively uniform 
effect. 

By similar reasoning, at most a factor two can be gained through parallelization of our 
implementation by assigning terms in the Hardy- Ramanuj an- Rademacher sum to separate 
threads. Further speedup on a multicore system requires parallelized versions of lower level 
routines, such as the exponential function or bignum multiplication. (A simplistic multithreaded 
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version of the FLINT partition function was tested for n up to 10 18 , giving nearly a twofold 
speedup on two cores, but failed when computing 10 19 for reasons yet to be determined.) 
Fortunately, it is likely to be more interesting in practice to be able to evaluate p{n) for a 
range of large values than just for a single value, and this task naturally parallelizes well. 



5. Multi-evaluation and congruence generation 

One of the central problems concerning the partition function is the distribution of values 
of p(n) mod m. In 2000, Ono |On00j proved that for every prime m > 5, there exist infinitely 
many congruences of the type 

p(Ak + B) = mod rn (5.1) 

where A, B are fixed and k ranges over all nonnegative integers. Ono's proof is nonconstructive, 
but Weaver |We01| subsequently gave an algorithm for finding congruences of this type when 
m G {13, 17, 19, 23, 29, 31}, and used the algorithm to compute 76,065 explicit congruences. 

Weaver's congruences are specified by a tuple (m, £, e) where £ is a prime and e G { — 1,0, 1}, 
where we unify the notation by writing (m,£, 0) in place of Weaver's (m,£). Such a tuple 
corresponds to a family of congruences of the form (|5.1|) with coefficients 

A = m£ 4 -^ (5.2) 
B = ^-^ + 1 + m^\'% (5.3) 

where a is the unique solution of m£ 3 ~^a = — 1 mod 24 with 1 < a < 24, and where < S < £ 
is any solution of 

(245 ^ -a mod £ if er = 

| (245 + a \£)=e ife = ±l. ^ " ' 

The free choice of 8 gives £ — 1 distinct congruences for a given tuple (m,£,s) if e = 0, and 
(£ — l)/2 congruences if e = ±1. 

Weaver's test for congruence, described by Theorems 7 and 8 in [WeOl], essentially 
amounts to a single evaluation of p(n) at a special point n. Namely, for given m,£, we 
compute the smallest solutions of S rn = 24 _1 mod m, r rn = —m mod 24, and check whether 
p(mr m (£ 2 -l)/24 + S m ) is congruent mod m to one of three values corresponding to the 
parameter e G { — 1,0,1}. We give a compact statement of this procedure as Algorithm [5] To 
find new congruences, we simply perform a brute force search over a set of candidate primes 
£, calling Algorithm [5] repeatedly. 



n 


Decimal 


expansion 


Number of digits 


Terms 


Error 


10 12 


6129000962 


. . 6867626906 


1,113,996 


264,526 


2 x 10~ 7 


10 13 


5714414687 


. .4630811575 


3,522,791 


787,010 


3 x 10~ 8 


10 14 


2750960597 


. . 5564896497 


11,140,072 


2,350,465 


-1 x 10- 8 


10 15 


1365537729 


. . 3764670692 


35,228,031 


7,043,140 


-3 x 10~ 9 


10 16 


9129131390 


. . 3100706231 


111,400,846 


21,166,305 


-9 x 10- 10 


10 17 


8291300791 


. . 3197824756 


352,280,442 


63,775,038 


5 x 10- 10 


10 18 


1478700310 


. . 1701612189 


1,114,008,610 


192,605,341 


4 x 10- 10 


10 19 


5646928403 


. . 3674631046 


3,522,804,578 


582,909,398 


4 x 10- 11 



Table 2. Large values of p(n). The table also lists the number of terms N in the 
Hardy-Ramanujan-Rademacher formula used by FLINT (theoretically bounding the error by 
0.25) and the difference between the floating-point sum and the rounded integer. 
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Algorithm 5 Weaver's congruence test 

Input: A pair of prime numbers 13 < ra < 31 and £ > 5, ra ^ £ 
Output: (m,£,s) defining a congruence, and Not-a-congruence otherwise 

8 rn ^— 24 _1 mod m {Reduced to < S m < ra} 

r m <- (-ra) mod 24 {Reduced to < ra < 24} 

V <r- ^ 

x <— p(S m ) {We have x ^ mod m} 

/ <- (3 | ^) ((-l)V m | {Jacobi symbols} 
t^y + fxt' 1 

if £ = cj mod m where cj G { — 1, 0, 1} then 

return (ra, cj (3(-l) v | £)) 
else 

return Not-a-congruence 
end if 



5.1. Comparison of algorithms for vector computation 

In addition to the Hardy-Ramanujan-Rademacher formula, the author has added code 
to FLINT for computing the vector of values p(0),p(l), . . .p{n) over Z and Z/mZ. The 
code is straightforward, simply calling the default FLINT routines for power series inversion 
over the respective coefficient rings, which in both cases invokes Newton iteration and FFT 
multiplication via Kronecker segmentation. 

A timing comparison between the various methods for vector computation is shown in 
Table [3l The power series method is clearly the best choice for computing all values up to 
n modulo a fixed prime, having a complexity of 0(n 1+ °^). For computing the full integer 
values, the power series and HRR methods both have complexity 0(n 3 / 2+ °^^), with the power 
series method expectedly winning. 

Ignoring logarithmic factors, we can expect the HRR formula to be better than the power 
series for multi-evaluation of p(n) up to some bound n when n/c values are needed. The factor 
c ~ 10 in the FLINT implementation is a remarkable improvement over c ~ 1000 attainable 
with previous implementations of the partition function. For evaluation mod ra, the HRR 
formula is competitive when (^(n 1 / 2 ) values are needed; in this case, the constant is highly 
sensitive to ra. 

For the sparse subset of C^n 1 / 2 ) terms searched with Weaver's algorithm, the HRR formula 
has the same complexity as the modular power series method, but as seen in Table [3] runs more 



n 


Series (Z/13Z) 


Series (Z) 


HRR (all) 


HRR (sparse) 


10 4 


0.01 s 


0.1 s 


1.4 s 


0.001 s 


10 5 


0.13 s 


4.1 s 


41 s 


0.008 s 


10 6 


1.4 s 


183 s 


1430 s 


0.08 s 


10 7 


14 s 






0.7 s 


10 8 


173 s 






8 s 


10 9 


2507 s 






85 s 



Table 3. Comparison of time needed to compute multiple values of p(n) up to the given 
bound, using power series inversion and the Hardy-Ramanujan-Rademacher formula. The 
rightmost column gives the time when only computing the subset of terms that are searched 
with Weaver's algorithm in the ra = 13 case. 
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than an order of magnitude faster. On top of this, it has the advantage of parallelizing trivially, 
being resumable from any point, and requiring very little memory (the power series evaluation 
mod m = 13 up to n = 10 9 required over 40 GiB memory, compared to a few megabytes with 
the HRR formula). Euler's method is, of course, also resumable from an arbitrary point, but 
this requires computing and storing all previous values. 

We mention that the authors of |CDJPS07] use a parallel version of the recursive Euler 
method. This is not as efficient as power series inversion, but allows the computation to be 
split across multiple processors more easily. 

5.2. Results 

Weaver gives 167 tuples, or 76,065 congruences, containing all £ up to approximately 1,000- 
3,000 (depending on m). This table was generated by computing all values of p(n) with n < 
7.5 x 10 6 using the recursive version of Euler's pentagonal theorem. Computing Weaver's table 
from scratch with FLINT, evaluating only the necessary n, takes just a few seconds. We are 
also able to numerically verify instances of all entries in Weaver's table for small k. 

As a more substantial exercise, we extend Weaver's table by determing all £ up to 10 6 for 
each prime m. Statistics are listed in Table HJ The computation was performed by assigning 
subsets of the search space to separate processes, running on between 40 and 48 active cores 
for a period of four days, evaluating p(n) at 6(7r(10 6 ) — 3) = 470, 970 distinct n ranging up to 
2 x 10 13 . 

We find a total of 70,359 tuples, corresponding to slightly more than 2.2 x 10 10 
new congruences. To pick an arbitrary, concrete example, one "small" new congruence is 
(13, 3797, -1) with S = 2588, giving 

p(711647853449/c + 485138482133) = mod 13 

which we easily evaluate for all k < 100, providing a sanity check on the identity as well as the 
partition function implementation. As a larger example, (29, 999959, 0) with S = 999958 gives 

p(28995244292486005245947069/c + 28995221336976431135321047) = mod 29 

which, despite our efforts, presently is out of reach for direct evaluation. 

Complete tables of (£, e) for each m are available at: 
http: //www. rise . jku. at/people/f johanss/partitions/ 
http : //sage .math. Washington. edu/home/f redrik/part it ions/ 



m 


(m,£,Q) 


(m,£,+l) 


(m,£,-l) 


Congruences 


CPU 


Max n 


13 


6,189 


6,000 


6,132 


5,857,728,831 


448 h 


5.9 x 10 12 


17 


4,611 


4,611 


4,615 


4,443,031,844 


391 h 


4.9 x 10 12 


19 


4,114 


4,153 


4,152 


3,966,125,921 


370 h 


3.9 x 10 12 


23 


3,354 


3,342 


3,461 


3,241,703,585 


125 h 


9.5 x 10 11 


29 


2,680 


2,777 


2,734 


2,629,279,740 


1,155 h 


2.2 x 10 13 


31 


2,428 


2,484 


2,522 


2,336,738,093 


972 h 


2.1 x 10 13 


All 


23,376 


23,367 


23,616 


22,474,608,014 


3,461 h 





Table 4. The number of tuples of the given type with £ < 10 6 , the total number of 
congruences defined by these tuples, the total CPU time, and the approximate bound up to 

which p(n) was evaluated. 



Page 18 of H9] 



FREDRIK JOHANSSON 



6. Discussion 



Two obvious improvements to our implementation would be to develop a rigorous, and 
perhaps faster, version of Algorithm [3] for computing cos(p7r/q) to high precision, and to develop 
fast multithreaded implementations of transcendental functions to allow computing p(n) for 
much larger n. Curiously, a particularly simple AGM-type iteration is known for exp(7r) (see 
[BB03]), and it is tempting to speculate whether a similar algorithm can be constructed for 
exp(7rV24n _ allowing faster evaluation of the first term. 

Some performance could also be gained with faster low-precision transcendental functions 
(up to a few thousand bits) and by using a better bound than ([1.8)) for the truncation error. 

The algorithms described in this paper can be adapted to evaluation of other HRR-type 
series, such as the number of partitions into distinct parts 



Using asymptotically fast methods for numerical evaluation of hypergeometric functions, it 
should be possible to retain quasi-optimality. 

Finally, it remains an open problem whether there is a fast way to compute the isolated 
value p(n) using purely algebraic methods. We mention the interesting recent work by Bruinier 
and Ono |BQ11| . which perhaps could lead to such an algorithm. 



The author thanks Silviu Radu for suggesting the application of extending Weaver's table 
of congruences, and for explaining Weaver's algorithm in detail. The author also thanks the 
anonymous referee for various suggestions, and Jonathan Bober for pointing out that Erdos' 
theorem about quadratic nonresidues gives a rigorous complexity bound without assuming the 
Extended Riemann Hypothesis. 

Finally, William Hart gave valuable feedback on various issues, and generously provided 
access to the computer hardware used for the large-scale computations reported in this paper. 
The hardware was funded by Hart's EPSRC Grant EP/G004870/1 (Algorithms in Algebraic 
Number Theory) and hosted at the University of Warwick. 
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